﻿using LxBasic;

namespace LxCore
{
    public class SIMPLEMass : MassFluxBase
    {
        public SIMPLEFields SimpleGroup;
        public override void CalFaceMassFlux(ThreadInfo threadInfo)
        {
            {
                (int i1, int endp1) = Zone.Inner.FaceIdx.Split(threadInfo);
                var den = SimpleGroup.Density.GetFaceVarArr(Zone.Inner);
                var vel = SimpleGroup.Velocity.GetInnCellVelVars();
                var InnDc = SimpleGroup.Pressure.InnDc;
                var GradP = SimpleGroup.Pressure.Grad;
                var P = SimpleGroup.Pressure.InnX;
                for (; i1 < endp1; i1++)
                {
                    ref FvInnFace face = ref InnFaces[i1];
                    Vector3d v_O = vel.GetCellVar(face.OSideCell);
                    Vector3d v_N = vel.GetCellVar(face.NSideCell);
                    Vector3d vf = (1 - face.gOwner) * v_O + face.gOwner * v_N;
                    {
                        var Dc_O = InnDc[face.OSideCell];
                        var Dc_N = InnDc[face.NSideCell];
                        var GradP_O = GradP[face.OSideCell];
                        var GradP_N = GradP[face.NSideCell];
                        var Df_ = (1 - face.gOwner) * Dc_O + face.gOwner * Dc_N;
                        var GradPf_ = (1 - face.gOwner) * GradP_O + face.gOwner * GradP_N;
                        var GradPf = GradPf_ + (P[face.NSideCell] - P[face.OSideCell] - Vector3d.DotProduct(GradPf_, face.OToN)) / (face.OToNLength * face.OToNLength) * face.OToN;
                        var DcGradPf_ = (1 - face.gOwner) * Matrix3dMath.Product(Dc_O, GradP_O) + face.gOwner * Matrix3dMath.Product(Dc_N, GradP_N);
                        vf += DcGradPf_ - Matrix3dMath.Product(Df_, GradPf);
                    }
                    double des = den.GetFaceVar(i1);
                    InnMassFlux[i1] = des * Vector3d.DotProduct(vf, face.Normal);
                }
            }
            var mbs = SimpleGroup.SimpleBounds;
            for (int i1 = 0; i1 < mbs.Length; i1++)
            {
                switch (mbs[i1])
                {
                    case SIMPLEFields.NoSlipWall noslip:
                        setWallFlux0(Zone.Bounds[i1], noslip, threadInfo);
                        break;
                    case SIMPLEFields.Symmetry sym:
                        setWallFlux0(Zone.Bounds[i1], sym, threadInfo);
                        break;
                    case SIMPLEFields.VelocityInlet velInlet:
                    case SIMPLEFields.PressureInOut p:
                        calFaceMassFlux_VelBound(Zone.Bounds[i1], threadInfo);
                        break;
                }
            }
            void setWallFlux0(FvBound b, SIMPLEFields.ZeroMassFlux mb, ThreadInfo threadInfo)
            {
                (int i1, int endp1) = b.FaceIdx.Split(threadInfo);
                Array.Clear(BouMassFlux, i1, endp1 - i1);
            }
            void calFaceMassFlux_VelBound(FvBound bound, ThreadInfo threadInfo)
            {
                (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
                var den = SimpleGroup.Density.GetFaceVarArr(bound);
                var vel = SimpleGroup.Velocity.GetBouCellVelVars();
                for (; i1 < endp1; i1++)
                {
                    ref var face = ref BouFaces[i1];
                    Vector3d vf = vel.GetCellVar(i1);//不进行修正，因为压力本身是根据内部网格计算出来的，不会出现棋盘型问题，且这样保证流入的质量一定等于速度乘以面积
                    BouMassFlux[i1] = den.GetFaceVar(i1) * Vector3d.DotProduct(vf, face.Normal);
                }
            }
        }
    }

}
